####################################################################
# Replication code for
#   Hirose, Imai, and Lyall, "Can civilian attitudes predict insurgent violence?"
#
# Code to produce IRT estimates of civilian support for ISAF and Taliban 
#
# Author: Kentaro Hirose
####################################################################



#############################################################################
chain <- 1 ### use four MCMC chains (i.e., chain = 1, 2, 3, 4) 
#############################################################################
library(endorse)
#############################################################################
load("survey.RData")
load("village.covariate.RData")
#############################################################################
index <- 1:length(unique(survey$PPLID))
data.village <- data.frame(index)
data.village$PPLID <- unique(survey$PPLID)
############################################################################
Y <- list(Q1 = c("Q5.1", "Q5.1A", "Q5.1B"),
          Q2 = c("Q5.2", "Q5.2A", "Q5.2B"),
          Q3 = c("Q5.3", "Q5.3A", "Q5.3B"),
          Q4 = c("Q5.4", "Q5.4A", "Q5.4B"))
############################################################################
formula.indiv <- formula(~ age +
                         educ +
                         madrassa +
                         pro.taliban +
                         pro.taliban.na + 
                         income +
                         income.na +
                         violent.exp.taliban +
                         violent.exp.taliban.na +
                         violent.exp.ISAF +
                         violent.exp.ISAF.na + 
                         encounter.taliban +
                         encounter.taliban.na +
                         encounter.ISAF +
                         encounter.ISAF.na +
                         frequency.ISAF)

formula.village <- formula(~ 1)
############################################################################
MCMC <- 200000
burn <- 180000
thin <- 40
set.seed(chain)
nu0.omega2 = 5
s0.omega2 = 2
nu0.phi2 = 5
s0.phi2 = 2
nu0.psi2 = 5
s0.psi2 = 2
s.start = rnorm(1)
x.start = rnorm(1)
############################################################################
endorse.out <- endorse(Y = Y,
                       na.strings = 99, 
                       data = survey,
                       data.village = data.village,
                       formula.indiv = formula.indiv,
                       formula.village = formula.village,
                       village = "PPLID",
                       identical.lambda = TRUE,
                       covariates = TRUE,
                       hierarchical = TRUE,
                       s.out = TRUE, 
                       x.sd = FALSE,
                       MCMC = MCMC,
                       burn = burn,
                       thin = thin,
                       nu0.omega2 = nu0.omega2,
                       s0.omega2 = s0.omega2,
                       nu0.phi2 = nu0.phi2, 
                       s0.phi2 = s0.phi2,
                       nu0.psi2 = nu0.psi2,
                       s0.psi2 = s0.psi2,
                       s.start = s.start,
                       x.start = x.start
                       )
############################################################################
x <- endorse.out$x
ip <- x
x <- as.matrix(apply(x, 1, sd))
s <- endorse.out$s/matrix(x, nrow(x), ncol(endorse.out$s))
############################################################################
nsim <- nrow(x)
S <- matrix(NA,nsim,nrow(survey))
for (i in 1:nsim){
  si <- s[i,]
  si <- matrix(si,4,)
  S[i,] <- apply(si,2,mean)
}
############################################################################
COND <- village.covariate$surveyed == 1 & !is.na(village.covariate$surveyed) & village.covariate$villageid %in% unique(survey$PPLID)
data.village <- village.covariate[COND,]
J <- nrow(data.village)
ideal.point <- support.isaf <- support.taliban <- matrix(NA, nsim, J)
verbose <- 1
for (j in 1:J){
  COND <- survey$endorsement.ISAF == 1 & survey$PPLID == data.village$villageid[j]
  A <- as.matrix(S[, COND])
  support.isaf[, j] <- apply(A, 1, mean)
  
  COND <- survey$endorsement.taliban == 1 & survey$PPLID == data.village$villageid[j]
  A <- as.matrix(S[, COND])
  support.taliban[, j] <- apply(A, 1, mean)

  COND <- survey$PPLID == data.village$villageid[j]
  A <- as.matrix(ip[, COND])
  ideal.point[, j] <- apply(A, 1, mean)
}
############################################################################
if(chain == 1){ ### MCMC results from chain 1
  write.csv(support.isaf, "support.isaf.1.csv")
  write.csv(support.taliban, "support.taliban.1.csv")
  write.csv(ideal.point, "ideal.point.1.csv")
}
if(chain == 2){ ### MCMC results from chain 2
  write.csv(support.isaf, "support.isaf.2.csv")
  write.csv(support.taliban, "support.taliban.2.csv")
  write.csv(ideal.point, "ideal.point.2.csv")
}
if(chain == 3){ ### MCMC results from chain 3
  write.csv(support.isaf, "support.isaf.3.csv")
  write.csv(support.taliban, "support.taliban.3.csv")
  write.csv(ideal.point, "ideal.point.3.csv")
}
if(chain == 4){ ### MCMC results from chain 4
  write.csv(support.isaf, "support.isaf.4.csv")
  write.csv(support.taliban, "support.taliban.4.csv")
  write.csv(ideal.point, "ideal.point.4.csv")
}
############################################################################
